function rie = computeRieImplicit(i,elements,model)
e = model.getElement(i-1);
f = e.getFace();
ng = e.getNumGaussPoint();
nn = f.countModes();
dim = 2;
rie = zeros(dim*nn,1);
for j=1:ng
    xi = e.getGaussPoint(j-1);
    w = e.getGaussWeight(j-1);
    jac = f.jacXAt(xi(1),xi(2));
    B = zeros(3,dim*nn);
    B = e.computeB(xi(1),xi(2),f,B);
    sigma = elements(i).stressImplicit(:,j);
    rie = rie + w*B'*sigma*abs(det(jac));
end
